% This code: assign a measure of road infrastructure to actual road network data

function roads = assign_weights_to_roads_LAC( roads )

% number of lanes
road_lanes_mat = cell2mat( {roads.lanes} );       % number of lanes
road_lanes_mat(isnan(road_lanes_mat)) = ones(sum(isnan(road_lanes_mat)), 1);

% intended use  % =0: local; =1: national
% classify highway, primary as national and all others as local
road_use_mat = cell2mat( {roads.primary} );
road_use_missing =  cell2mat( {roads.missing_pr} );     % missing lanes 
road_use_mat(road_use_missing==1) = zeros(sum(isnan(road_use_mat)), 1);
 
% paved? lanes % =0: unpaved; =1: paved
road_paved_mat = cell2mat( {roads.paved} );
road_paved_missing = isnan(road_paved_mat);
road_paved_mat(isnan(road_paved_mat)) = ones(sum(isnan(road_paved_mat)), 1);

% distance
road_length_mat = deg2km( cell2mat( {roads.Shape_Leng} ) );

% cheapest-path coefficients
chi_lane = 0.2;
chi_use = 1.07;
chi_paved = 1.35;

% weight
shape_weight = road_length_mat.*...
                ( road_lanes_mat.^( -chi_lane ) ).*...
                 chi_use.^( 1-road_use_mat ).*...
                 chi_paved.^( 1-road_paved_mat );   % weight of entire shape

%% assign a distance in km and a weight to each edge

for j=1:length(roads)
    
    [ j length(roads) ]
    
    roads(j).totdist = road_length_mat(j); % deg2km( roads(j).SHAPE_Leng ); % total distance, in KM
    if isempty(roads(j).totdist)
        roads(j).totdist=0; 
    end
    
    % length of each link in km
    Ntemp = length(roads(j).X)-1;
    x0 = roads(j).X(1:Ntemp-1);
    y0 = roads(j).Y(1:Ntemp-1);
    x1 = roads(j).X(2:Ntemp);
    y1 = roads(j).Y(2:Ntemp);  
    if roads(j).totdist>0
        roads(j).distances = distance( y0,x0,y1,x1 )./sum( distance( y0,x0,y1,x1 ) ) * roads(j).totdist;       
    elseif roads(j).totdist==0
        roads(j).distances = distance( y0,x0,y1,x1 );
    end
    
    % link features
    roads(j).lanes = road_lanes_mat(j);
    roads(j).totlanes = roads(j).lanes*roads(j).totdist;    
    roads(j).use = road_use_mat(j);
    roads(j).median = 1;
    roads(j).paved = road_paved_mat(j);
    
    % missing values
    roads(j).missing_lanes = road_lanes_mat(j)==1;
    roads(j).missing_use = road_use_missing(j);
    roads(j).missing_paved = road_paved_missing(j);
    
    % segment weight
    roads(j).weights = distance( y0,x0,y1,x1 )./sum( distance( y0,x0,y1,x1 ) ) * shape_weight(j);
    
end
